#!/usr/bin/perl -w
use strict;
use Getopt::Long;
use Bio::Perl;
use Bio::DB::Fasta;
use Bio::SeqIO;
use List::Util qw(min max);

my $flank = 5000;
my $minmatch = 70;
my $name = 'TE';
my $dir  = ".";
my $debug = 0;
my $Min_evalue = 1e-30;
GetOptions(
	   'f|flank:s'   => \$flank,
	   'n|name:s'    => \$name,
	   'o|d|dir:s'   => \$dir,
	   'v|verbose!'  => \$debug,
	   );
mkdir($dir) unless -d $dir;
# order of argument is
# Chromosome FASTA database file (genome)
# Proteins searched with (pepdb)
# Table of BLAST results in -mformat 3 from WU-BLAST
my ($genome,$pepdb,$table) = @ARGV;

my $pepdbh = Bio::DB::Fasta->new($pepdb);
my $dbh = Bio::DB::Fasta->new($genome);

open(my $fh => $table) || die $!;
my (%groups,%segments);
my $ct = 0;
# parse BLAST (WU-BLAST mformat=3 and -links), key part is the -links option
while(<$fh>) {
    next if /^\#/;
    my ($qid,$sid, $evalue, $N, $Sprime,$S,$alnlen,$nident,$npos,$nmism,$pcident,$pcpos,
	$qgaps,$qgaplen,$sgaps,$sgaplen,$qframe,$qstart,$qend,
	$sframe,$sstart,$send, $links) = split;
   # next if $pcpos < $minmatch;
    my $num = 1;
    if ($links =~ /\((\d+)\)/) {
	($num) = $1;
    }
    $links =~ s/[\(\)]//g;
    $groups{$sid}->{$qid}->{$links} = $evalue;
    $segments{$sid}->{$qid}->[$num] = [$sstart => $send,$pcpos];    
    $ct++;
}
$ct = 0;
# walk through all the chromosomes with hits
for my $chrom ( keys %groups ) {
    for my $q ( keys %{$groups{$chrom}} ) { # walk through all the proteins with hits
	    for my $grp ( keys %{$groups{$chrom}->{$q}} ) { # walk through each of these logical groups (this is gotten by the -links option in WU-BLAST)
	    next if $groups{$chrom}->{$q}->{$grp} > $Min_evalue;
	    my ($min,$max,$low_pcident);
	    for my $num ( split(/-/,$grp) ) {
		my $hsp = $segments{$chrom}->{$q}->[$num];		
		if( ! defined $hsp->[0] ) {
		    warn("no value for $grp:$chrom:$q\n");
		}
		$min = ! defined $min ? $hsp->[0] : min($hsp->[0],$min);
		$max = ! defined $max ? $hsp->[1] : max($hsp->[1],$max);
		if( $hsp->[2] < $minmatch ) {
		    $low_pcident = 1;
		}
	    }
	    if( ! defined $min ) {
		warn("no value for min for grp:$grp in $chrom-$q\n");
		next;
	    }
	    my ($left,$right) = ($min - $flank, 
				 $max   + $flank);
	    my $chromlen = $dbh->length($chrom);
	    $left = 1 if $left < 1;
	    $right = $chromlen if $right > $chromlen;

	    my $segment = $dbh->seq($chrom,$left => $right);	    
	    my $oname = "$dir/$name\_$ct";
	    my $out = Bio::SeqIO->new(-format => 'fasta',
				      -file   => ">$oname.fa");
	    $out->write_seq(Bio::Seq->new(-id => "$name\_$ct",
					  -desc => sprintf("%s:%d..%d target=%d..%d domain=%s",$chrom,$left,$right,$min,$max,$q),
					  -seq => $segment));
	    my $qid = $q;
	    $qid =~ s/\/(\d+)\-(\d+)$//;
	    if( ! -f "$dir/$qid" && ! -z "$dir/$qid" ) {
		   $out = Bio::SeqIO->new(-format => 'fasta',-file => ">$dir/$qid")->write_seq(Bio::Seq->new(-id => $qid,
													  -seq => $pepdbh->seq($q)));
	    }
      # Do a realignment of the protein to genome to handle this better than the BLAST alignment. This is supposed to basically be replacement of the PHI heuristic with 
      # a full dynamic programming alignment
	    open(my $ofh => "exonerate -m protein2genome -q $dir/$qid -t $oname.fa --refine full --ryo \"##startCDS\n>%ti %tab..%tae (%tcl nt)\n%tcs\n##endCDS\n\" |") || die $!;
	    my $cds;
	    my $keep = 0;
	    while(<$ofh> ) {
		if(/^\#\#startCDS/ ) {
		   $keep = 1;
	       } elsif( /^\#\#endCDS/ ) {
		   $keep = 0;
	       } elsif( $keep ) {		       
		   $cds .= $_;
	       }		
	    }
	    my $cds_seq = Bio::SeqIO->new(-fh => IO::String->new($cds),-format => 'fasta')->next_seq;
	    $out = Bio::SeqIO->new(-format =>'fasta',-file => ">$oname.pep");
	    $out->write_seq($cds_seq->translate);
	    $ct++;
	    last if $debug;
	}
	last if $debug;
    }
    last if $debug;
}

# collapses runs of numbers into groups
sub collapse_nums {
# This is probably not the slickest connectivity algorithm, but will do for now.
    my @a = @_;
    my ($from, $to, $i, @ca, $consec);
    
    $consec = 0;
    for($i=0; $i < @a; $i++) {
        not $from and do{ $from = $a[$i]; next; };
    # pass repeated positions (gap inserts)
    next if $a[$i] == $a[$i-1];
        if($a[$i] == $a[$i-1]+1) {
            $to = $a[$i];
            $consec++;
        } else {
            if($consec == 1) { $from .= ",$to"; }
            else { $from .= $consec>1 ? "\-$to" : ""; }
            push @ca, split(',', $from);
            $from =  $a[$i];
            $consec = 0;
            $to = undef;
        }
    }
    if(defined $to) {
        if($consec == 1) { $from .= ",$to"; }
        else { $from .= $consec>1 ? "\-$to" : ""; }
    }
    push @ca, split(',', $from) if $from;

    @ca;
}

